🎓 Homework #1 🎓

---

👤 Assigner:

Malkova Ksenia

✍️ Contacts:

• kemalkova@edu.hse.ru

• @ksu_marshmallow

🔥 Deadline 🔥:

• 11.02.2025

---

import os
import warnings
import subprocess
from pathlib import Path

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA

import seaborn as sns
import plotly.express as px
import matplotlib.pyplot as plt

sns.set(style="whitegrid", context="notebook", palette="pastel")
warnings.simplefilter('ignore')

DATA_FOLDER = Path('../data/HW1')
SCRIPT_FOLDER = Path('../scripts')
SEED = 42

🧬 Task 1

Для набора из трех позиций и четырех сэмплов, сделайте PCA. Спроектируйте на главные компоненты, мотивировав чем-то выбор количества компонент.

Матрица:

$$X = \begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} $$

Решение:

1. Центрируем матрицу $$ Z = X - \bar{X} $$

• Средние по столбцам: $\bar{X_1} = \bar{X_2} = \bar{X_3} = 0.25$
• Стандартизировать не будем, все и так в одном масштабе

$$ Z = \begin{pmatrix} -0.25 & 0.75 & -0.25 \\ -0.25 & -0.25 & 0.75 \\ 0.75 & -0.25 & -0.25 \\ -0.25 & -0.25 & -0.25 \\ \end{pmatrix} $$

---

2. Строим ковариационную матрицу

Ковариация (несмещенная оценка): $$ \text{cov}(X, Y) = \frac{1}{n-1} \sum_i^n (X - \bar{X}) \cdot (Y - \bar{Y}) $$

Ковариационная матрица - матрица, у которой $(i,j)$-элемент является ковариацией признаков $(X_i, X_j)$. Когда $X_i = X_j$ ковариация - это дисперсия $X_j$ $(\mathbb{D}X_j)$

• По диагонали - дисперсии признаков
• В остальных ячейках - ковариации соответствующих признаков
• Т.к. ковариация симметричная, то и матрица ковариации симметрична

$$ \text{COV\_MATRIX} = \begin{pmatrix} \text{cov}(X_1, X_1) & \text{cov}(X_1, X_2) & \text{cov}(X_1, X_3) \\ \text{cov}(X_2, X_1) & \text{cov}(X_2, X_2) & \text{cov}(X_2, X_3) \\ \text{cov}(X_3, X_1) & \text{cov}(X_3, X_2) & \text{cov}(X_3, X_3) \\ \end{pmatrix} $$

Наша ковариационная матрица: $$ \mathbb{D}X_1 = \mathbb{D}X_2 = \mathbb{D}X_3 = \frac{(1 + 1 + 1 + 9)}{3 \cdot 4^2} = \frac{1}{4} $$

$$ \text{cov}(X_1, X_2) = \text{cov}(X_1, X_3) = \text{cov}(X_2, X_3) = \frac{-3 + 1 - 3 + 1}{3 \cdot 4^2} = -\frac{1}{12} $$

$$ \text{COV} = \begin{pmatrix} \frac{1}{4} & -\frac{1}{12} & -\frac{1}{12} \\ -\frac{1}{12} & \frac{1}{4} & -\frac{1}{12} \\ -\frac{1}{12} & -\frac{1}{12} & \frac{1}{4} \\ \end{pmatrix} $$

---

3. Ищем собственные числа и вектора в ковариационной матрице.

Собственные значения находят через уравнение (3), после чего с известными с.з. находят собственные вектора через уравнение (2) $$ \begin{gather} A \cdot v = \lambda \cdot v \\ v \cdot (A - \lambda I) = 0 \\ det(A - \lambda I) = 0 \end{gather} $$ где $A$ - исходная матрица, $v$ - собственный вектор, $\lambda$ - собственное значение, $I$ - identity matrix

• Каждому собственному вектору соответствует свое собственное значение (пары собственные вектор-значение)
• Собственные значения отражают "важность" компоненты - какую долю дисперсии изначальных признаков она несет (больше дисперсия $\to$ больше информации содержит данный признак. Т.е. чем больше собственное значение, тем важнее компонента)
• Каждый собственный вектор показывает направление, вдоль которого наблюдается соответствующая дисперсия. Это и есть компоненты в новом пространстве признаков!
• Компоненты в новом пространстве независимы друг от друга

$$ \begin{gather} det(A - \lambda I) = 0 \\ \begin{vmatrix} \frac{1}{4} - \lambda && -\frac{1}{12} && -\frac{1}{12} \\ -\frac{1}{12} && \frac{1}{4} - \lambda && -\frac{1}{12} \\ -\frac{1}{12} && -\frac{1}{12} && \frac{1}{4} - \lambda \\ \end{vmatrix} = 0 \\ \end{gather} $$

Определитель матрицы 3*3: $$ \begin{gather} (\frac{1}{4} - \lambda)^3 - 2 \cdot \frac{1}{12^3} - 3 \cdot \frac{1}{12^2} \cdot (\frac{1}{4} - \lambda) = 0 \\ \frac{(3-12\lambda)^3 - 2 - 3 \cdot (3 - 12 \lambda)}{12^3} = 0 \\ (3-12\lambda)^3 - 3 \cdot (3 - 12 \lambda) - 2 = 0 \\ \end{gather} $$

Пусть $t = 3 - 12 \lambda$ (т.е. $\lambda = \frac{3 - t}{12}$): $$ t^3 - 3 t - 2 = 0$$

Первый корень (перебором): $t_1 = -1 \to \lambda_1 = \frac{1}{3}$

$$ \begin{gather} (t + 1) (t^2 - t - 2) = 0 \\ t^2 - t - 2 = 0 \\ D = 1 + 8 = 9 \\ t_2 = \frac{1-3}{2} = -1 \to \lambda_2 = \frac{1}{3} \\ t_3 = \frac{1+3}{2} = 2 \to \lambda_3 = \frac{1}{12} \end{gather} $$

Собственные числа: $\frac{1}{12}$, $\frac{1}{3}$

Для каждого с.ч. найдем свой собственный вектор: $$v (A - \lambda I) = 0$$

1. $\lambda_1 = \frac{1}{12}$:

$$ \begin{gather} \begin{cases} \frac{2}{12} x_1 -\frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 + \frac{2}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 -\frac{1}{12} x_2 + \frac{2}{12} x_3 = 0 \\ \end{cases} && \to v_1 = (1, 1, 1) \end{gather} $$

2. $\lambda_2 = \frac{1}{3}$:

$$ \begin{gather} \begin{cases} -\frac{1}{12} x_1 -\frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 - \frac{1}{12} x_2 -\frac{1}{12} x_3 = 0 \\ -\frac{1}{12} x_1 -\frac{1}{12} x_2 - \frac{1}{12} x_3 = 0 \\ \end{cases} && \to v_2 = (-1, 1, 0), && v_3 = (-1, 0, 1) \end{gather} $$

Пары собственные числа-собственные значения: $$ \begin{gather} \lambda_1 = 1/12 && v_1 = (1, 1, 1) \\ \lambda_2 = 1/3 && v_2 = (-1, 1, 0) \\ \lambda_3 = 1/3 && v_3 = (-1, 0, 1) \end{gather} $$

---

4. Проекция на главные компоненты

Главные компоненты - топ-K компонент, несущих бОльшую часть информацию (имеющих бОльшие собственные значения). С одной стороны, мы хотим снизить размерность данных, а с другой - сохранить больше информации.

Возьмем топ-$K$ компонент, которые кумулятивно объясняют ~90% вариации исходных данных. (На практике же, когда признаков очень много, число $K$ определяют "методом локтя"=elbow method)

Сначала отсортируем собственные вектора в соответствие с их собственными значениями:

$$\Lambda = (\lambda_2, \lambda_3, \lambda_1) = \big( \frac{1}{3}, \frac{1}{3}, \frac{1}{12} \big)$$

$$ V = (v_2, v_3, v_1) = \begin{pmatrix} -1 & -1 & 1 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \\ \end{pmatrix} $$

И найдем долю объясненной дисперсии $d$, которую несет каждая компонента - по сути каждое собственное значение делим на сумму всех собственных значений.

$$ \begin{gather} d_2 = d_3 = \frac{1/3}{1/3+1/3+1/12} = \frac{4}{9} \\ d_1 = \frac{1/12}{1/3+1/3+1/12} = \frac{1}{9} \end{gather} $$

Чтобы понять, до какого числа компонент сокращать - должны посчитать долю дисперсии кумулятивно:

$$ \begin{gather} d_{PC_1} = 4/9 \\ d_{PC_2} = 8/9 \\ d_{PC_3} = 1 \end{gather} $$

Двумя главными компонентами можем объяснить 88.88 $\approx$ 90% дисперсии исходных данных. Поэтому третью компоненту скипаем, оставляем первые две. И далее делаем проекцию исходных данных на главные компоненты. Проецируя данные, мы переводим их из исходной системы координат (где переменные могут быть коррелированы) в новую, где каждая компонента независима от других.

Математически - ищем dot-product между исходными (центрированными) данными $Z$ и матрицей базисных векторов (топ-$K$ компонент) $V_K$: $$X' = Z \cdot V_K$$ $$ \begin{gather} X' = \begin{pmatrix} -0.25 & 0.75 & -0.25 \\ -0.25 & -0.25 & 0.75 \\ 0.75 & -0.25 & -0.25 \\ -0.25 & -0.25 & -0.25 \\ \end{pmatrix} \cdot \begin{pmatrix} -1 && -1 \\ 1 && 0 \\ 0 && 1 \\ \end{pmatrix} \\ X ' = \begin{pmatrix} 1 && 0 \\ 0 && 1 \\ -1 && -1 \\ 0 && 0 \end{pmatrix} \end{gather} $$

---

И код прилагается, но с.в. и трансформированная матрица будут различаться

# Initialize matrix
X = np.array([
    [0, 1, 0],
    [0, 0, 1],
    [1, 0, 0],
    [0, 0, 0]
])

# 1. Centralize matrix
means = X.mean(axis=0)
stds = np.sqrt(X.var(axis=0))

print(f'Means: {means.round(2)} \nStandard deviations: {stds.round(2)}')

Z = (X - means) # / stds
Z

Means: [0.25 0.25 0.25] 
Standard deviations: [0.43 0.43 0.43]

array([[-0.25,  0.75, -0.25],
       [-0.25, -0.25,  0.75],
       [ 0.75, -0.25, -0.25],
       [-0.25, -0.25, -0.25]])

# 2. Compute covariance matrix
cov_matrix = np.cov(Z, rowvar=False)
cov_matrix

array([[ 0.25      , -0.08333333, -0.08333333],
       [-0.08333333,  0.25      , -0.08333333],
       [-0.08333333, -0.08333333,  0.25      ]])

# 3. Compute eigenvalues and eigenvectors
eigen_values, eigen_vectors = np.linalg.eig(cov_matrix)

print(f"Eigen-values: {eigen_values} \n\nEigen-vectors: \n{eigen_vectors.round(2)}")

Eigen-values: [0.33333333 0.08333333 0.33333333] 

Eigen-vectors: 
[[ 0.82 -0.58 -0.22]
 [-0.41 -0.58 -0.57]
 [-0.41 -0.58  0.79]]

# 4. Sort eigenvalues and eigenvectors in descending order of eigenvalues
order = np.argsort(eigen_values)[::-1]
eigen_vectors = eigen_vectors[:, order]
eigen_values = eigen_values[order]

# Find explained variance ratio to choose number of components
explained_variance_ratio_ = eigen_values / np.sum(eigen_values)
explained_variance_ratio_.cumsum().round(2)

array([0.44, 0.89, 1.  ])

# 5. Project data into new space 
n_components = 2
X_transformed = X.dot(eigen_vectors[:, :n_components])

X_transformed.round(4)

array([[-0.5709, -0.4082],
       [ 0.791 , -0.4082],
       [-0.2201,  0.8165],
       [ 0.    ,  0.    ]])

---

🧬 Task 2

Дано:

1. Файлы с вариантами chr21.pca.txt, chr22.pca.txt
2. Текстовый файл со списком сэплов IBS.YRI.MEX.txt
3. Файл со списком сэплов и их принадлежностью к популяция IBS.YRI.MEX.info.txt
4. Ноутбук pop.model.l3.ipynb c PCA по 22ой хромосоме для популяций IBS, YRI, MEX.

def read_data(file, f_names, f_info):
    with open(f_names, 'r') as f:
        sample_names = [line.strip() for line in f]
    
    with open(f_info, 'r') as f:
        info_data = [line.strip().split('\t') for line in f]
    
    info_dict = {row[0]: (row[1], row[3]) for row in info_data}

    populations = [info_dict[name][1] for name in sample_names]

    with open(file, 'r') as f:
        genotypes = [
            [
                0 if allele == '0|0' else 2 if allele == '1|1' else 1
                for allele in line.strip().split('\t')[3].split()
            ]
            for line in f
        ]

    genotype_matrix = np.array(genotypes).T

    genotype_df = pd.DataFrame(genotype_matrix)
    meta_df = pd.DataFrame({
        'Pop': populations,
        'ID': sample_names
    })

    result_df = pd.concat([genotype_df, meta_df], axis=1)
    return result_df

Какой процент дисперсии описывает первая главная компонента для 22ой хромосомы?

gt22 = read_data(file=DATA_FOLDER / 'chr22.pca.txt',
                 f_names=DATA_FOLDER / 'IBS.YRI.MEX.txt',
                 f_info=DATA_FOLDER / 'IBS.YRI.MEX.info.txt'
                 )

print(gt22.shape)
gt22.head()

(278, 18271)

	0	1	2	3	4	5	6	7	8	9	...	18261	18262	18263	18264	18265	18266	18267	18268	Pop	ID
0	0	2	1	0	1	0	2	1	0	2	...	0	1	0	0	0	1	0	0	IBS	HG01503
1	1	1	0	1	1	1	1	1	0	2	...	0	2	0	1	0	1	0	0	IBS	HG01510
2	0	1	0	0	1	0	1	0	0	0	...	0	1	0	0	0	0	0	0	IBS	HG01515
3	0	2	0	0	2	0	2	0	1	1	...	0	0	0	0	0	1	0	0	IBS	HG01522
4	1	1	1	1	1	1	1	1	0	0	...	0	0	0	0	1	0	0	0	IBS	HG01527

5 rows × 18271 columns

X22 = gt22.drop(['Pop', 'ID'], axis=1)

pca22 = PCA(n_components=3)
X22_transformed = pca22.fit_transform(X22)

explained_variance_ratio_22 = pca22.explained_variance_ratio_
print(f'Первая компонента объясняет ~{100 * explained_variance_ratio_22[0]:.2f}% вариации данных')

Первая компонента объясняет ~15.64% вариации данных

plt.figure(figsize=(8, 6), dpi=80)

sns.barplot(x = [f'PC{x}' for x in range(1, 4)], y = explained_variance_ratio_22,
            color='skyblue')

plt.xlabel('Principal components')
plt.ylabel('Explained variance ratio');

plt.figure(figsize=(8, 6), dpi=120)

sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'].values)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 22 chr genotype data");

---

Постройте PCA по 21ой хромосоме. Есть ли существенные отличия от визуализации по 22ой хромосоме?

gt21 = read_data(file=DATA_FOLDER / 'chr21.pca.txt',
                 f_names=DATA_FOLDER / 'IBS.YRI.MEX.txt',
                 f_info=DATA_FOLDER /  'IBS.YRI.MEX.info.txt'
                 )

print(gt21.shape)
gt21.head()

(278, 18860)

	0	1	2	3	4	5	6	7	8	9	...	18850	18851	18852	18853	18854	18855	18856	18857	Pop	ID
0	1	1	2	1	0	2	1	2	1	2	...	1	1	1	1	2	1	0	2	IBS	HG01503
1	0	1	2	1	1	1	1	1	2	1	...	2	0	1	0	1	1	1	1	IBS	HG01510
2	0	1	1	0	0	0	1	1	0	1	...	2	0	1	0	2	2	0	2	IBS	HG01515
3	1	1	1	1	0	1	1	1	1	1	...	1	1	2	0	2	2	0	1	IBS	HG01522
4	1	0	2	0	1	1	1	2	1	2	...	2	0	1	0	1	1	1	1	IBS	HG01527

5 rows × 18860 columns

X21 = gt21.drop(['Pop', 'ID'], axis=1)

pca21 = PCA(n_components=3)
X21_transformed = pca21.fit_transform(X21)

plt.figure(figsize=(8, 6), dpi=120)

sns.scatterplot(x=X21_transformed[:, 0], y=X21_transformed[:, 1], hue=gt21['Pop'].values)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 21 chr genotype data");

Визуально отличий нет, YRI сгруппировались в кучу где-то отдельно (сильно отличаются от двух других популяций). MXL и IBS близки по первой компоненте, но разделяются второй. Также между ними есть небольшое пересечение

Спроектируйте наблюдения, построенных по 22ой хромосоме, на главные вектора, построенных по 21ой хромосоме, сравните с рисунком, где данные 22ой хромосомы спроектированы на вектора, полученных по 22ой хромосоме.

У 21 и 22 хромосом разное количество снипов. Не вижу другого выхода, кроме как обрезать по минимальному числу, иначе отображение не получится.

print(f'Features num: \nchr21 - {gt21.shape[1]}, chr22 - {gt22.shape[1]}')

Features num: 
chr21 - 18860, chr22 - 18271

num_features = min(X21.shape[1], X22.shape[1])

X21_filtered = X21.iloc[:, :num_features]
X22_filtered = X22.iloc[:, :num_features]

assert X21_filtered.shape == X22_filtered.shape

Переделываем PCA 21-й хромосомы, так как фичи изменились

pca21_new = PCA(n_components=3)
pca21_new.fit(X21_filtered)

PCA(n_components=3)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Проецируем наблюдения с 22-й хромосомы на вектора 21-й хромосомы

X22_transformed_from_21 = X22_filtered.dot(pca21_new.components_.T).to_numpy()
X22_transformed_from_21.shape

(278, 3)

Рисуем ✨картиночки✨

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,6), dpi=120, sharex=True, sharey=True)

sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_transformed_from_21[:, 0], y=X22_transformed_from_21[:, 1], hue=gt22['Pop'], ax=ax2)

ax1.set_xlabel('PC1')
ax2.set_xlabel('PC1')

ax1.set_ylabel('PC2')

ax1.set_title('PCA on 22 chr')
ax2.set_title('PCA on 22 chr with 21 chr principal components');

Красота 💅

Объясняется тем, что PCA для 21-й хромосомы выделяет направления наибольшей дисперсии именно для снипов этой хромосомы. Главные компоненты 21-й хромосомы не отражают вариацию 22-й хромосомы, так как ее "ковариационная структура" выглядит иначе. Соответствия между важными позициями 21-й и 22-й хромосом, скорее всего нет (разве что на уровне случайности). Допустим, что в 21-й хромосоме важнее (более вариабельные) позиции {X}, а в 22-й - {Y}, и эти позиции не перекрываются. Главные компоненты для 21-й хромосомы будут переводить данные в пространство, где будет сохраняться вариабельность именно этих позиций. Если в 22-й хромосоме на этих позициях низкая вариабельность, или ее нет, то в новом пространстве ее и не станет. При этом для вариабельных позиций 22-й хромосомы при переходе в новое пространство исходня вариабельность затеряется, так как в матрице главных компонент 21-хромосомы этой структуры не было. В результате наблюдения 22-й хромосомы оказываются сжатыми в точку, так как по выбранным направлениям различий практически нет.

---

Объедините варианты 21ой и 22ой хромосом. Постройте PCA. Сравните результат с предыдущими.

X = pd.concat((X21, X22), axis=1)
X.columns = np.concatenate((X21.columns.values, 
                            X22.columns.values + X21.columns.max()))
print(X.shape)
X.head()

(278, 37127)

	0	1	2	3	4	5	6	7	8	9	...	37116	37117	37118	37119	37120	37121	37122	37123	37124	37125
0	1	1	2	1	0	2	1	2	1	2	...	1	1	0	1	0	0	0	1	0	0
1	0	1	2	1	1	1	1	1	2	1	...	2	0	0	2	0	1	0	1	0	0
2	0	1	1	0	0	0	1	1	0	1	...	1	0	0	1	0	0	0	0	0	0
3	1	1	1	1	0	1	1	1	1	1	...	0	0	0	0	0	0	0	1	0	0
4	1	0	2	0	1	1	1	2	1	2	...	0	0	0	0	0	0	1	0	0	0

5 rows × 37127 columns

pca = PCA(n_components=3)
X_transformed = pca.fit_transform(X)

X_transformed.shape

(278, 3)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6), dpi=120, sharex=True, sharey=True)

sns.scatterplot(x=X21_transformed[:, 0], y=X21_transformed[:, 1], hue=gt21['Pop'], ax=ax1)
sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax2)
sns.scatterplot(x=X_transformed[:, 0], y=X_transformed[:, 1], hue=gt22['Pop'], ax=ax3)

for ax, name in zip((ax1, ax2, ax3), ('chr21', 'chr22', 'chr21 + chr22')):
    ax.set_xlabel('PC1')
    ax.set_ylabel('PC2')

    ax.set_title(f'PCA on {name}')

plt.show()

Особо ничего не меняется. Только увеличилось "внутригрупповое" разнообразие внутри MXL - они стали более отдален друг от друга и от IBS (хотя некоторые наблюдения все равно пересекаются с этой популяцией)

---

Нарисуйте трехмерную проекцию на три главных компоненты для 22 хромосомы.

def plot_3d(x: np.ndarray, y: np.ndarray, z: np.ndarray, 
            populations: np.ndarray, 
            title: str = None, save_name: str = None
            ):

    fig = px.scatter_3d(
        x=x, y=y, z=z,
        color=populations,
        title=title,
        opacity=0.8
    )
    
    fig.update_traces(marker=dict(size=4))
    
    fig.update_layout(
        scene=dict(
            xaxis_title="PC1",
            yaxis_title="PC2",
            zaxis_title="PC3"
        ),
        legend=dict(
            title="Population",
            x=1.05, y=1,
            bgcolor="rgba(255, 255, 255, 0.8)"
        ),
        margin=dict(l=0, r=0, b=0, t=30)
    )
    
    fig.show()

    if save_name is not None:
        fig.write_html(f"{save_name}.html")

plot_3d(x=X22_transformed[:, 0], 
        y=X22_transformed[:, 1], 
        z=X22_transformed[:, 2], 
        populations=gt22['Pop'].values, 
        title='PCA on chr22',
        save_name='../imgs/HW1_task2_chr22'
        )

---

Случайным образом уменьшите количество маркеров на 50 и на 80 процентов. Визуализируйте и проанализируйте.

Допустим это касается только 22-й хромосомы

# 1. Отбираем 50 и 20% случайных столбцов
X22_20 = X22.sample(frac=0.2, axis=1, random_state=SEED)
X22_50 = X22.sample(frac=0.5, axis=1, random_state=SEED)

X22.shape[1], X22_20.shape[1], X22_50.shape[1]

(18269, 3654, 9134)

# 2. Делаем pca
pca_22_20 = PCA(n_components=3)
pca_22_50 = PCA(n_components=3)

X22_20_transformed = pca_22_20.fit_transform(X22_20)
X22_50_transformed = pca_22_50.fit_transform(X22_50)

# 3. Визуализируем
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6), dpi=120, sharex=True, sharey=True)

sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_50_transformed[:, 0], y=X22_50_transformed[:, 1], hue=gt22['Pop'], ax=ax2)
sns.scatterplot(x=X22_20_transformed[:, 0], y=X22_20_transformed[:, 1], hue=gt22['Pop'], ax=ax3)

for ax, name in zip((ax1, ax2, ax3), ('chr22', 'chr22 (50% markers)', 'chr22 (20% markers)')):
    ax.set_xlabel('PC1')
    ax.set_ylabel('PC2')

    ax.set_title(f'PCA on {name}')

plt.show()

Как при уменьшении до 50%, так и до 20% количества маркеров, заметно разделение популяций, хотя оно становится более "сжатым". Чем больше изначальных признаков мы удаляем, тем больше теряем информации о дисперсии - можем утратить важные направления, вдоль которых наблюдается вариация. Тем не менее, оказалось, что (особенно при 50%) это не так критично для описания различий между тремя популяциями, ведь они все равно различимы (хотя без раскраски это не было бы так очевидно, особенно для IBS и YRI при сокращении до 20% маркеров).

Структура данных сохраняется, вероятно, потому, что популяционные различия отражаются не отдельными маркерами а группами связанных маркеров. То есть многие маркеры несут избыточную информацию о популяционных различиях, поэтому удаление их части не критично, пока остаются маркеры, отражающие те же самые различия. Однако при сильном уменьшении количества маркеров вероятность того, что потеряются и такие группы, возрастает, что ведет к уменьшению "разрешающей способности" PCA

---

Уравняйте количество в каждой популяции. Изменится ли от этого PCA?

gt22.Pop.value_counts()

Pop
YRI    108
IBS    106
MXL     64
Name: count, dtype: int64

Сделаем клиппинг по наименее распространенной популяции MXL. Скорее всего результат не изменится - не сильно выражена доминантность какой-то популяции (-ий).

pop_count = gt22.Pop.value_counts().min().item()

gt22_balanced = gt22.groupby('Pop').sample(pop_count, random_state=SEED)
X22_balanced = gt22_balanced.drop(['Pop', 'ID'], axis=1)

pca_22_balanced = PCA(n_components=3)
X22_balanced_transformed = pca_22_balanced.fit_transform(X22_balanced)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,6), dpi=120, sharex=True, sharey=True)

sns.scatterplot(x=X22_transformed[:, 0], y=X22_transformed[:, 1], hue=gt22['Pop'], ax=ax1)
sns.scatterplot(x=X22_balanced_transformed[:, 0], y=X22_balanced_transformed[:, 1], hue=gt22_balanced['Pop'], ax=ax2)

ax1.set_xlabel('PC1')
ax2.set_xlabel('PC1')

ax1.set_ylabel('PC2')

ax1.set_title('PCA on 22 chr')
ax2.set_title('PCA on 22 chr (balanced populations)');

Ничего существенно не изменилось, разве что цвета точек поменялись, что связано с изменившимся порядком в таблице.

Это указывает на то, что в изначальных данных преобладают межпопуляционные различия. То есть различия между популяциями выражены значительно сильнее, чем вариация внутри каждой популяции. Даже после удаления какого-то количества образцов из двух популяций (почти половина!) межгрупповая вариация сохраняется. Поскольку межгрупповая вариация доминирует над внутригрупповой, PCA ее захватывает и отображает в первых компонентах. Скорее всего, изменились только последние главные компоненты, которые отвечают за малые внутрипопуляционные различия.

---

🧬Task 3

1. Выберите 4 ваших любимых популяции из проекта 1000GP Phase 3.

   • Список популяций и сэмплов можно найти тут
   • Скачать .vcf файлы chr22 можно здесь

2. Сделайте два текстовых файла: один со списком сэмплов, расположенных в столбик; второй состоит из двух столбцов Sample, Population.

3. Скачайте и запустите скрипт pca.sh. Пропишите в нем путь до переменной VCF0 и путь до plink. Запустив его, вы получите текстовый файл, в котором данные отфильтрованы.

4. Нарисуйте проекцию на две главные компоненты. Сопоставьте с тем, что известно про эти популяции из истории и географии.

BASE_URL = 'https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/'

VCF_URL = BASE_URL + 'ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz'
ANNOTATION_URL = BASE_URL + 'integrated_call_samples_v3.20200731.ALL.ped'

VCF_PATH = DATA_FOLDER /  '1KGP_phase3_22chr.vcf.gz'
SAMPLE_PATH = DATA_FOLDER /  '1KGP_phase3_22chr_samples'
ANNOTATION_PATH = DATA_FOLDER /  '1KGP_phase3_22chr_annotation'

1. Скачиваем .vcf 22-й хромосомы и ее индекс

subprocess.run(f'wget {VCF_URL} -O {VCF_PATH} -q', shell=True)
subprocess.run(f'wget {VCF_URL}.tbi -O {VCF_PATH}.tbi -q', shell=True)

CompletedProcess(args='wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz.tbi -O ../data/HW1/1KGP_phase3_22chr.vcf.gz.tbi -q', returncode=0)

2. Скачиваем список sample_id, отбираем из них только те, которые соответствуют выбранным популяциям
   • GBR - British in England and Scotland
   • CHB - Han Chinese in Beijing, China
   • LWK - Luhya in Webuye, Kenya
   • JPT - Japanese in Tokyo, Japan

Решила усложнить себе жизнь

Еще, чтобы воспользоваться функцией чтения файлов, я делала файл не sample-популяция, а sample-gender-sample-population (я не поняла, что в 3 колонке должно быть, но она и не используется)

samples_script ="""
#!/bin/bash

DATA_FOLDER=$HOME/PopGen_Models/data/HW1

ANNOTATION_URL="https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20200731.ALL.ped"
ANNOTATION_PATH=$DATA_FOLDER/1KGP_phase3_22chr_annotation
SAMPLE_PATH=$DATA_FOLDER/1KGP_phase3_22chr_samples

# Get all samples
bcftools query -l ${DATA_FOLDER}/1KGP_phase3_22chr.vcf.gz > all_samples

# Install annotation
wget ${ANNOTATION_URL} -q -O tmp.csv

# Query populations -> 'SAMPLE-POPULATION' file
awk '$7 == "GBR" || $7 == "CHB" || $7 == "LWK" || $7 == "JPT" {print $2 "\t" $5 "\t" $2 "\t" $7}' tmp.csv > annotation_samples

# Intersect with .vcf samples
join -1 1 -2 1 -t $'\t' <(sort all_samples) <(sort annotation_samples) > ${ANNOTATION_PATH}
echo "Annotation saved in ${ANNOTATION_PATH}"

# 'SAMPLE' file
awk '{print $1}' ${ANNOTATION_PATH} > ${SAMPLE_PATH}
echo "Samples saved in ${SAMPLE_PATH}"

rm tmp.csv annotation_samples all_samples
"""

SAMPLE_EXEC_PATH = SCRIPT_FOLDER / 'get_samples.sh'

with open(SAMPLE_EXEC_PATH, 'w') as f:
    print(samples_script, file=f)

subprocess.run(f'chmod +x {SAMPLE_EXEC_PATH}', shell=True)
subprocess.run(f'bash {SAMPLE_EXEC_PATH}', shell=True)

Annotation saved in /home/ksu_marshmallow/PopGen_Models/data/HW1/1KGP_phase3_22chr_annotation
Samples saved in /home/ksu_marshmallow/PopGen_Models/data/HW1/1KGP_phase3_22chr_samples

CompletedProcess(args='bash ../scripts/get_samples.sh', returncode=0)

3. Запускаем скрипт pca.sh. Я в нем только изменила пути, чтобы ко мне в папку с данными все скачивалось + изменила название выходного файла, чтобы не перезаписывался файл из второй таски.

PCA_EXEC_PATH = SCRIPT_FOLDER / 'pca.sh'

subprocess.run(f'chmod +x {PCA_EXEC_PATH}', shell=True)
subprocess.run(f'{PCA_EXEC_PATH} 22 {SAMPLE_PATH}', shell=True)

PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.log.
Options in effect:
  --allow-extra-chr
  --double-id
  --indep-pairwise 50 10 0.1
  --out /home/ksu_marshmallow/PopGen_Models/data/HW1/test2
  --set-missing-var-ids @:#
  --vcf /home/ksu_marshmallow/PopGen_Models/data/HW1/test.vcf.gz

15374 MB RAM detected; reserving 7687 MB for main workspace.
--vcf: /home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.bed +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.bim +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test2-temporary.fam written.
1054487 variants loaded from .bim file.
1054487 missing IDs set.
397 people (0 males, 0 females, 397 ambiguous) loaded from .fam.
Ambiguous sex IDs written to
/home/ksu_marshmallow/PopGen_Models/data/HW1/test2.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 397 founders and 0 nonfounders present.
Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done.
1054487 variants and 397 people pass filters and QC.
Note: No phenotypes present.
Pruned 770911 variants from chromosome 22, leaving 283576.
Pruning complete.  770911 of 1054487 variants removed.
Marker lists written to
/home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.in and
/home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.out .
PLINK v1.90b6.21 64-bit (19 Oct 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.log.
Options in effect:
  --allow-extra-chr
  --double-id
  --extract /home/ksu_marshmallow/PopGen_Models/data/HW1/test2.prune.in
  --maf 0.1
  --make-bed
  --out /home/ksu_marshmallow/PopGen_Models/data/HW1/test3
  --set-missing-var-ids @:#
  --vcf /home/ksu_marshmallow/PopGen_Models/data/HW1/test.vcf.gz

15374 MB RAM detected; reserving 7687 MB for main workspace.
--vcf: /home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.bed +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.bim +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test3-temporary.fam written.
1054487 variants loaded from .bim file.
1054487 missing IDs set.
397 people (0 males, 0 females, 397 ambiguous) loaded from .fam.
Ambiguous sex IDs written to
/home/ksu_marshmallow/PopGen_Models/data/HW1/test3.nosex .
--extract: 283576 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 397 founders and 0 nonfounders present.
Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done.
266392 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
17184 variants and 397 people pass filters and QC.
Note: No phenotypes present.
--make-bed to /home/ksu_marshmallow/PopGen_Models/data/HW1/test3.bed +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test3.bim +
/home/ksu_marshmallow/PopGen_Models/data/HW1/test3.fam ... 101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899done.

CompletedProcess(args='../scripts/pca.sh 22 ../data/HW1/1KGP_phase3_22chr_samples', returncode=0)

4. Делаем PCA

gt22_task3 = read_data(file=DATA_FOLDER / 'task3_chr22.pca.txt',
                       f_names=SAMPLE_PATH,
                       f_info=ANNOTATION_PATH
                       )

print(gt22_task3.shape)
gt22_task3.head()

(397, 17186)

	0	1	2	3	4	5	6	7	8	9	...	17176	17177	17178	17179	17180	17181	17182	17183	Pop	ID
0	1	1	0	1	2	1	1	0	0	0	...	1	1	0	1	0	0	1	0	GBR	HG00096
1	0	1	0	0	0	0	1	0	0	1	...	1	1	1	1	1	0	1	0	GBR	HG00097
2	0	1	0	0	0	0	1	0	0	1	...	0	1	0	0	0	1	0	0	GBR	HG00099
3	0	2	1	0	1	0	2	0	1	2	...	1	0	0	0	0	2	0	0	GBR	HG00100
4	0	2	1	0	1	0	2	1	1	2	...	1	1	0	1	0	2	0	0	GBR	HG00101

5 rows × 17186 columns

# Распределение популяций
gt22_task3.Pop.value_counts()

Pop
JPT    104
CHB    103
LWK     99
GBR     91
Name: count, dtype: int64

X22_task3 = gt22_task3.drop(['Pop', 'ID'], axis=1)

pca22_task3 = PCA(n_components=3)
X22_task3_transformed = pca22_task3.fit_transform(X22_task3)

plt.figure(figsize=(8, 6), dpi=120)

sns.scatterplot(x=X22_task3_transformed[:, 0], y=X22_task3_transformed[:, 1], hue=gt22_task3['Pop'].values)

plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title("PCA on 22 chr genotype data");

Давным-давно, в самой древней Африке, жили люди. Это было время, когда все люди на Земле были как одна большая семья. Они жили рядом друг с другом и знали только один дом - африканские равнины.

Однажды эта семья решила отправиться в большое путешествие. Кто-то остался дома, в теплых и солнечных краях (это были предки наших кенийцев, LWK), а другие захотели увидеть мир. И вот они начали свой долгий путь...

Некоторые из них пошли на северо-запад, через горы и пустыни. Дорога была нелегкой, но они нашли новые земли - холодные и суровые. Эти смельчаки стали предками европейцев, а потом и британцев (GBR). В этих новых местах им пришлось учиться жить в холоде, одеваться теплее и находить новую еду. Их тела и гены начали меняться, чтобы приспособиться к новому (дивному) миру.

А другие отправились совсем в другую сторону - на восток. Они шли через Индию, через горы и реки, пока не добрались до великой Азии. Там они тоже разделились: одни остались в Китае (CHB), а другие через какое-то время переплыли море и обосновались на островах Японии (JPT). Хотя они теперь жили далеко друг от друга, их предки когда-то были одной семьей, поэтому их гены все еще очень похожи.

Те, кто остался в Африке (LWK), продолжали жить так же, как их предки. Они знали каждый уголок своей родной земли, приспособились к ее жаре и дождям. Их гены оставались особенными, потому что они никогда не покидали свой дом.

Спустя тысячи лет, когда ученые начали изучать гены людей со всего мира, они заметили удивительную вещь. Гены тех, кто живет в Европе (британцы), совсем не похожи на гены тех, кто живет в Африке (кенийцы). А гены китайцев и японцев, хотя они и разные, все еще напоминают друг друга - ведь их предки когда-то были одной семьей.

И вот, когда мудрые ученые построили свою волшебную карту PCA, она показала три большие семьи: африканцы (LWK) остались в своем древнем доме, европейцы (GBR) ушли далеко на север, а азиаты (CHB и JPT) отправились на восток. Первые две главные дороги на карте (PC1 и PC2) рассказали о большом разделении между континентами. А китайцы с японцами, хотя и живут отдельно уже давно, все еще стоят близко — их история расхождения спрятана где-то на более тонких тропинках, может быть, на третьей компоненте.

Так генетика помогла нам прочитать древнюю историю человечества, записанную в наших клетках 🌍✨

plot_3d(x=X22_task3_transformed[:, 0], 
        y=X22_task3_transformed[:, 1], 
        z=X22_task3_transformed[:, 2], 
        populations=gt22_task3['Pop'].values, 
        title='PCA on chr22',
        save_name='../imgs/HW1_task3_chr22'
        )

Действительно! Третья компонента PCA уже может разделять японцев и китайцев. Хоть тоже не идеально и они пересекаются, но хотя бы не накладываются друг на друга, как только с первыми двумя векторами.